Thermalization and ergodicity in one-dimensional many-body open quantum systems 
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Using an approach based on the time-dependent density-matrix renormalization group method, we 
study thermalization in spin chains locally coupled to an external bath. Our results provide evidence 
that quantum chaotic systems do thermalize, that is, they exhibit relaxation to an invariant ergodic 
state which, in the bulk, is well approximated by the grand canonical state. Moreover, the resulting 
ergodic state in the bulk does not depend on the details of the baths. On the other hand, for 
integrable systems we found that the invariant state in general depends on the bath and is different 
from the grand canonical state. 
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The emergence of canonical ensembles in quantum sta- 
tistical mechanics from first principles is one of the key 
remaining old questions of theoretical physics. Even the 
definition of the temperature at the nano-scale poses a 
challenge ;lj]. Namely, the main question is how to "de- 
rive" the canonical distribution? It has been realized 
that the canonical distribution is in a way "typical" : 
provided the overall system describing the environment 
plus a central system is in a generic pure state, the re- 
duced state of the central system is with high probabil- 
ity canonical [2|. However, how precisely the canonical 
distribution arises from dynamical laws, without a pri- 
ori statistical assumptions, is still unclear. Motivation in 
the study of this fundamental aspect of nonequilibrium 
physics also comes from some recent experiments with ul- 
tracold bosonic gases, where absence of thermalization in 
closed, integrable, strongly correlated quantum systems 
has been observed Q. 

For closed many-body systems, integrability is believed 
to play a crucial role in the relaxation to the Steady State 
(SS): the nonequilibrium dynamics of a chaotic system is 
expected to thermalize at the level of individual eigen- 
states H, as numerically observed in several physical 
models [5|. By contrast, for systems with non trivial 
integrals of motion, SSs usually carry memory of the ini- 
tial conditions and are not canonical: maximizing the 
entropy while keeping the values of constants of motion 
fixed results in a generalized Gibbs ensemble ^6]. Much 
less is known about the relaxation to the SS for open 
quantum systems this is what we are going to ad- 
dress in this paper. We provide numerical evidence that, 
analogously to closed systems, the occurrence of thermal- 
ization is strictly related to system's integrability, irre- 
spective of the fine details of the baths. In particular we 
show that locally coupling a quantum chaotic many-body 
system to an environment is enough for a SS of the cen- 



tral system to be very close, in the bulk, to the canonical 
or grand canonical state (GCS). On the contrary, if the 
system is integrable, the constants of motion in general 
prevent thermalization and the form of the SS sensitively 
depends on the bath coupling operators. We show that 
the numerical description of an open quantum system in 
terms of a Lindblad equation with local coupling to the 
reservoirs is in some sense a computationally efficient, 
minimal model of thermalization. Such result paves the 
way for future simulations of quantum transport in large 
many-body quantum systems. 

The time evolution for a generic state p of an open 
quantum system can be described, under certain approx- 
imations, by a Lindblad master equation Q: 

^P- |[p,H]+£bP, (1) 

where Ti is the Hamiltonian of the autonomous system, 
while the dissipation Cb = 7j2k {[^kP^Ll] + [Lk,pLl]j 

is parametrized by certain Lindblad operators Lk (here- 
after we set h = Ub = ^ and, unless noted otherwise, 
7 = 1). The derivation of Eq. ([T|) from first principles, 
i.e., from the Hamiltonian evolution of a system plus envi- 
ronment is rather tricky Q ; however it is the most general 
form of a completely positive, trace preserving, dynam- 
ical semi-group. Taking it for granted, we ask ourselves 
if, within this approximation, a finite many-body system 
can thermalize when coupled via some Lindblad opera- 
tors Lk acting only locally just on few degrees of freedom. 

To elucidate the role covered by chaoticity in the 
thermalization process, we consider prototype one- 
dimensional spin-l/2-chain models with nearest neighbor 
interactions: % — Ym=o {hi^i+i denoting the local 

energy density, and n being the chain length). As we 
shall see, the chosen models exhibit a crossover from in- 
tegrable to chaotic regime when a suitable parameter in 
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their Hamiltonians is varied. With the term "chaotic" we 
refer, as usual, to a system whose bulk energy spectrum 
of highly excited levels obeys a random matrix statis- 
tics Q; in particular, the level spacing statistics (LSS) 
p{s) is well approximated by the Wigner-Dyson distri- 
bution pwd(s) @j whereas in an integrable system LSS 
typically turns out to be Poissonian, pp{s). 

We assume local coupling to the reservoirs, i.e., the 
dissipator Cb acts only on the m n) leftmost (/) and 
rightmost (r) spins: jCb — J^s'^^buik'SijCQ. We construct 
Cb by generalizing the method discussed in Ref. ■ For 
this purpose, we first consider the GCS for the spin chain, 



pg{T, fi) = exp [^{n - fi E^)/r] 



(2) 



where = J2?=o '^i total magnetization [cr" 

{a = x, y,z) being the Pauli operators for the jth spin], 
T the temperature, ^ the "chemical potential", and 
Z = tr [exp (— — I]^)/T)] the partition function. 
Given a target temperature Ttarg and a chemical po- 
tential /itarg, the reduced m-spin target density matrix 
Ptkrg: ^ £ {^^}: is obtaiucd after tracing (Ttarg, Mtarg) 
over all but the m leftmost /rightmost spins. We finally 
require that p^^rg the unique eigenvector of £3 with 
eigenvalue 0, while all other eigenvalues are equal to 
— 1. Such a choice produces, in absence of % and for 
a given spectral norm of /^g, the fastest convergence to 

Ptarg 



In the presence of % we obtain, for up to 



n « 100 spins, the SS solution of Eq. ([T]) numerically 
by using a time-dependent Density Matrix Remormal- 
ization Group (tDMRG) method with a Matrix Product 
Operator (MPO) ansatz fli| . 

In the following we are interested in the asymptotic 
state reached, independently of initial conditions, after a 
long time, pss = hmt-i-oo p{t)- In all simulations we care- 
fully checked that the simulation time was long enough 
to reach convergence, which is exponential. Since Lind- 
blad operators act only locally and pg (T, p) is invariant 
for the unitary part of Eq. ([T]), pss cannot be equal to 
the GCS, unless it is also an eigenstate of the dissipator 
Cb- In other words, one can have pss = Pg{T, p) only if 
pg(r, p) = pj^jg ® pbulk ® pLrg: i-e-, if the GCS is sepa- 
rable with respect to the border m spins which are used 
in the coupling. Nevertheless for chaotic systems, as we 
shall see, sufficiently far from the boundaries the state is 
arbitrarily close to pg (T, p) , regardless of the entangle- 
ment with the coupled parts. 

Let us start our numerical investigations by consid- 
ering a spin-1/2 Ising chain in a tilted magnetic field, 
described by the energy density 



hi. 



)■ (3) 



Its only conserved quantity is the total energy, therefore 
the expected invariant state is the canonical one pg{T, 0). 
To check thermalization, we solved the master equation 
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FIG. 1: (Color online). One-spin observables for the SS (sym- 
bols) agree with the theoretical canonical ones (full curve) 
to less than 0.5% for the chaotic Ising model (main plot). 
For the integrable model (inset) a comparable agreement is 
observed by looking at energy density (since (cr^^j) = 0) 
vs. {(7^/2) ■ Squares are for n = 16 and uniform couplings; cir- 
cles (triangles) for n = 16 (n = 40) and couplings Ji switched 
on over a layer of thickness r = 4, with 7 = 0.2. Marks on 
theoretical curves show the temperature. 



for two different sets of parameters: (i) a transverse field 
6x 1, &z = 0, for which the model is integrable and ex- 
hibits a Poissonian LSS; (ii) a tilted field = 1, &z = 1, 
for which it is chaotic with a Wigner-Dyson LSS [ISj (if 
not specified, we take Ji — 1 and couple two border spins, 
m — 2). With the obtained pss, '^e evaluated expecta- 
tion values of several one- and two-spin observables in the 
bulk of the chain, and compared them to the theoretical 
ones as given by the canonical state pg(T, 0). 

In the main plot of Fig. [T] we show one-spin expecta- 
tion values (c^^j) — tr {pss f "72) fo^' the chaotic case: all 
numerical points fall on the curve given by theoretical ex- 
pectation values for a canonical state. The same happens 
in the integrable Ising model. Such irrelevance of inte- 
grability is a peculiarity of certain few-body observables, 
similarly to what observed in a different context of out- 
of-equilibrium dynamics in closed systems [l^ . Quite 
remarkably, we could not reach temperatures in the bulk 
below « 1.7 (see squares in Fig. [1]), even by using very 
small Ttarg ~ 0. The reason resides in the already men- 
tioned boundary effects due to entanglement between the 
boundary two spins and the bulk chain, which makes the 
cooling difficult. This must be contrasted with a zero at- 
tainable temperature in the case of separable states . 
For entangled states though, our results show that to 
lower the minimal attainable temperature one has to re- 
duce the effect of interaction at the boundaries which is 
responsible for entanglement. One way to do this is by 
switching on the interaction gently over a boundary layer 
of certain thickness r, Ji = sm{-^^) {Jn-2-1 = sin{-^^)), 
for Z = 0, . . . , r — 1, at the left (right) end and using a 
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FIG. 2: (Color online). Absolute differences in the expecta- 
tion value of a two-spin observable o-^/2'^n/2+i between the 
SS and the theoretical canonical state, in the case of chaotic 
(full symbols) and integrable Ising model (empty symbols). 
Dashed lines denote constant relative error. 
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FIG. 3: (Color onhne). SS (symbols) and GCS (full lines) 
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expectation values of ' (left panel), a 
(right panel) for the Heisenberg model with n — ti't) spins, 
Ttarg = 4, Qtarg = 2, m = 3, J; = 1. lu the left panel, "XX" 
and "XXZ" refer to two integrable cases without magnetic 
field (respectively at A = 0, 0.5), "stagg." to the chaotic case 
with A = 0.5 and period-3 staggered field with B — 2. The 
curves in the right panel are for the chaotic case only. The 
GCS ps(rmoas ~ 5.851, /imcas ~ —0.534) for the chaotic case 
is obtained by matching (/i3i+i,3i+2) and {crli^i) for which 
lines are not shown in the right panel. 



weaker coupling 7 (circles and triangles in Fig.[T]). 

To make comparison between pss and pg (T, 0) quanti- 
tative, we determined the "measured" temperature Tmcas 
to which Pss corresponds, which is in general different 
from Ttarg, due to boundary effects. Assuming that the 
SS is canonical in the bulk, one can extract Tmcas by 
comparing observables that uniquely set the tempera- 
ture. For Ising model (jSj, the energy density is sufficient, 
therefore we used the condition tr [ft,„/2-i.n/2 Pss] = 
tr [/i„/2-i,„/2 Pe(T'mcas,0)] to computc Tnicas- We then 
calculated theoretical expectation values of other observ- 
ables, through (Tineas 7 0); a comparison with the corre- 
sponding values for the reached SS may serve as an indi- 
cator of the quality of thermalization. In Fig. [5] we show 
differences between expectation values of afaf^-^^, com- 
puted with Pss and pc;(T,„cas, 0), for both chaotic and 
integrable Ising chains. A marked distinction between 
the two cases appears. First, in the chaotic model er- 
rors are much smaller than in the integrable one; second, 
switching J; gradually, which should decrease errors due 
to smaller boundary effects, in the integrable case even 
worsens the situation. The integrable Ising model there- 
fore does not relax to a canonical state in the bulk. Simi- 
lar results are obtained for other few-spin observables, as 
well as for the lowest moments of the energy distribution: 
we evaluated {[{'HQ — {'He))/5]P) {p = 2, 5 and He is the 
Hamiltonian of the 6 central spins) on the states pss and 
Ps(Tmcas,0). In a chain of n = 40 spins relative errors 
are never greater than 1% in the chaotic case, and are 
typically an order of magnitude larger in the integrable 
case. 

To corroborate the importance of system's integrability 
on the convergence to invariant statistical ensembles, we 
consider another prototype model of interacting spins: 
the Heisenberg XXZ chain in a magnetic field, described 



by the energy density 

(4) 

If the field is homogeneous the model is integrable and 
possesses, besides energy and magnetization, an infinite 
sequence of conserved quantities jl6|. On the other hand, 
integrability can be broken, e.g., simply by means of a 
period-3 staggered magnetic field, b^k = —B, b^k+i — 
—Bjl, bsk+2 — 0. In order to highlight the lack of ther- 
malization in the integrable regime B = 0, we target 
a non-Gibbsian state different from the GCS; namely, 
we use pnon-e(T, g) exp{-H/T + qQ^), with (34 = 

'^J '^i+i'^i+2'^J+3^ being a conserved charge for an open 
chain with A = and 6/ = [3]. The idea is that, 
using a 9targ 7^ 0, in the integrable regime the SS ex- 
hibits strong deviations from the GCS, corresponding to 
g = 0, while we expect chaotic dynamics to drive the bulk 
towards the GCS. Such expectation is confirmed by our 
numerical data. In Fig.|3]we show the spatial dependence 
of various observables for integrable, as well as for chaotic 
cases. In the integrable cases deviations from the GCS 
expectations are large, while they become very small for 
a chaotic system. Analogously to the Ising model, we 
checked this statement also for other few-spin observ- 
ables (we found that, in presence of chaos, the largest 
discrepancy among all the one- and two-spin observables 
amounts to 2 x 10"''^); layered interactions in the inte- 
grable model do not help in thermalizing the system. 

A further confirmation of the role of integrability 
comes from a direct analysis of the quality of thermaliza- 
tion after gradually switching on the perturbation that 
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FIG. 4: (Color online). Relative differences Ag^.^,/ = 
Ag<*'(B)/Ag<*'(0) in g'"** expectation values on the SS and 
the GCS evaluated in the bulk of the Heisenberg model 
with A = 0.5, as the staggering strength B is varied (full 
curves). Also shown is dependence of the rj function (dashed 
curves), characterizing the integrable-chaotic crossover. For 
both quantities the crossover takes place at smaller B with 
increasing n. Inset: two examples of LSS in the integrable 
{B = 0.01) and chaotic (B = 1) regimes; dashed lines denote 
Poissonian and Wigner-Dyson statistics fH, respectively. 



drives the crossover from integrability to chaos: the lon- 
gitudinal field bz in Eq. ^ or the staggering intensity B 
in Eq. As shown in Fig. |3] for the Heisenberg model, 
such crossover is conveniently detected by the parameter 
V = / -pwD(s)|ds// |pp(s) -pwD(s)|ds; 77 = 1 
and r] — correspond to Poissonian and Wigner-Dyson 
distributions, respectively. In the same figure we also plot 
deviations in (q*-"'-') evaluated on the SS and on the corre- 
sponding GCS: Aq^^'' = tr [(j|'*''(pss - pg(rmcas, /imoas))] 

as the strength B of the staggered magnetic field is in- 
creased. The progressive onset of chaos gradually im- 
proves the quality of thermalization, being Aq^'^^ a mono- 
tonic decreasing function of B. Moreover, the strength 
of the staggered field required to converge to the GC ex- 
pectation value drops with the system size. 

In conclusion, we have shown that, within the Lindblad 
equation formalism, coupling a one-dimensional quantum 
chaotic system locally to a bath results in a SS being 
equal to the invariant (grand) canonical state, far away 
from the coupled sites. In contrast, integrable systems 
do not thermalize and their SSs exhibit strong devia- 
tions from the (grand) canonical state, depending on the 
details of the coupling. The fact that for chaotic systems 
the SS does not depend on the details of the coupling, 
shows that very likely the same result would be obtained 
even for a harder-to-treat Hamiltonian evolution of a sys- 
tem plus environment or for higher dimensional systems. 
Our method should be applicable also to non-equilibrium 
situations. Indeed, by locally coupling a system to two 
or several baths at different values of temperature and 
chemical potentials, one should be able to efficiently con- 



trol local thermalization. Thus, our results might open 
significant new perspectives in the simulation of quan- 
tum transport in many-body quantum systems in contact 
with thermal and chemical baths. 
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